%this file is to calculate the 2-D correlation of u via FFT; it takes the 
%fluctuation fields as input

%Created by Ricardo Mejia-Alvarez.
%University of Illinois at Urbana-Champaign
%Urbana, IL. 06/16/2010

%modified:  06/22/2010

clear

%input data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
display('how many different data folders are you going to consider?');
NF = input(': ');
folData = cell(NF);
pathFiles = cell(NF);
NOR = zeros(NF);

for k = 1 : NF
   
    kstr = num2str(NF);
    message = strcat('Enter the path for the folder number ',kstr,' \n');
    folData{k} = input( message , 's');
    if ispc == 1
    folData{k} = strcat(folData{k},'\');
    else
        folData{k} = strcat(folData{k},'/');
    end
    
end

NameResult = input('Enter a name for your results: ','s');

display('Now, for each data folder, tell me what columns to use to calculate')
display('two-point correlations (use vector format).')

for k = 1 : NF
   
    kstr = num2str(NF);
    message = strcat('Columns for data folder number ',kstr,' \n');
    DataColumns{k} = input( message ); %#ok<*SAGROW>
    cols(k) = length(DataColumns{k}); 
end

cols = sum(cols);
disp('What is the extension of the data files?')
disp('1. dat')
disp('2. V3D')
NK = input(': ');

if NK == 1
    namekey = '*.dat';
else
    namekey = '*.V3D';
end

disp('Enter data transformation if necessary in matrix form.')
disp('For instance: ')
disp('                    X  Y  U1  U2  U3  ...  Un')
disp('Inversion         [-1  1  -1   1   1  ...   1 ]')
disp('Horizontal flip   [ 1  0   1   0   0  ...   0 ]')
disp('Vertical flip     [ 0  1   0   0   0  ...   0 ]')
disp('Transpose         [ 0  1   0   0   0  ...   0 ]')
disp('  ')
disp('Enter only the cells to modify.')
disp('E.g. T(1,1)=-1; T(1,3)=-1; T(2,1)=1; T(2,3)=1; T(3,2)=1; T(4,2)=1')
disp('   ')
disp('Type "return" and press the Return key to finish')
T = ones(1 , cols + 2);
T(2 : 4 , :) = zeros(3 , cols + 2);

keyboard

WRInt = input('Enter the wall range in vector format [starting_row  ending_row]');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if matlabpool('size') >= 1
    matlabpool close
end
matlabpool
sizePool = matlabpool('size');
dist = codistributor1d(3);

% %extent of reference range
WR = WRInt(end) - WRInt(1) + 1;

%creating the directory structure of the analysis
[~,folResults,folFluc] = newfolder(folData{1},'results','fluctuations'); 

for k = 1 : NF
    
    %generating a vector with the name of the files to process
    pathFiles{k} = filesVector(folData{k},namekey);

    %finding the number of instantaneous realizations
    NOR(k) = length(pathFiles{k});

end

testLength = isequal(NOR , repmat(NOR(1) , [1 , length(NOR)] ) );
if testLength == 0
    display('The number of files in the different folders is not the same.')
    display('I think you should take a look at this.')
    display('hit Return to continue')
    pause
    clear testLength
end

NOR = NOR(1);

%defining a minimum number of iterations (this will be redefined according
%to the size of the field)
NIT = 2;
nn = 1;
tic
while nn <= NIT
    mm = (nn - 1) * WR + 1;
    
    if nn <= NIT - 1
        WRInt = [mm , mm + WR - 1];
    else
        WRInt = [mm , NRPF];
        WR = NRPF - mm + 1;
    end
        
    tic
    for n = 1 : NOR

        p = 0;
        for k = 1 : NF
            filename = char(pathFiles{k}(n));
            [~,I,J,Dx,Dy,A] = matrix(filename);

            for m = 1 : length( DataColumns{k} )

                U(p + m) = { A{ DataColumns{k}(m) } };

            end
            p = length(U);

            [U] = Transformation(T(: , 3 : end) , U);

            if k == 1 && n == 1
                [A(1 : 2)] = Transformation(T(: , 1 : 2) , A(1 : 2));

                X = A{1};
                Y = A{2};

                if T(4 , 1) ~= 1
                    NCPF = I;
                    NRPF = J;
                else
                    NCPF = J;
                    NRPF = I;
                    Y = X;
                    Dx = Dy;
                end
                
                clear X S
            end

        end

        if n == 1

            %determine the length for fft transform
            Length = 2 * 2^(ceil(log(NCPF)/log(2)));
            L0 = Length - 2 * NCPF;

            %finding the horizontal and vertical "lag" coordinates
            Rx = Dx * [1 - NCPF : 0 , 1 : NCPF - 1] ;
            Ry = Y(: , 1)  ;
        end

        [SpectraCell , Sigma_YRefcell , Sigma_Ycell] = Spectra(U , WRInt , WR , NRPF , Length);

        if n ~= 1
            for m = 1 : length(SpectraCell)
                Spec_rho(m) = { Spec_rho{m} + SpectraCell{m} };  
            end

            for m = 1 : length(U)
                Sigma_YRef(m) = { Sigma_YRef{m} + Sigma_YRefcell{m} };  
                Sigma_Y(m) = { Sigma_Y{m} + Sigma_Ycell{m} };  
            end
        else

            Spec_rho = SpectraCell;
            Sigma_YRef = Sigma_YRefcell;
            Sigma_Y = Sigma_Ycell;
        end


    end
    toc

    disp('calculating correlations')
    tic

    for m = 1 : cols
        Sigma_YRef(m) = { Sigma_YRef{m} / ( NOR * Length ) };
        Sigma_Y(m) = { Sigma_Y{m} / ( NOR * Length ) };


        %Calculate RMS for y_ref
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        Sigma_YRef(m) = { gather( Sigma_YRef{m} ) };
        Sigma_YRef(m) = { real( ifft( Sigma_YRef{m} , Length , 2 ) ) };
        Sigma_YRef(m) = { sqrt( Sigma_YRef{m} ) };
%         Sigma_YRef(m) = { gather( Sigma_YRef{m} ) };
        Sigma_YRef(m) = { repmat( real( Sigma_YRef{m}(: , 1 , :) ) , [1 , 2 * NCPF - 1 , 1]) };
        Sigma_YRef(m) = { codistributed(Sigma_YRef{m} , dist) };
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


        %Calculate RMS for y
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        Sigma_Y(m) = { gather( Sigma_Y{m} ) };
        Sigma_Y(m) = { real( ifft( Sigma_Y{m} , Length , 2) ) };
        Sigma_Y(m) = { sqrt( Sigma_Y{m} ) };
%         Sigma_Y(m) = { gather( Sigma_Y{m} ) };
        Sigma_Y(m) = { repmat( real( Sigma_Y{m}(: , 1 , :) ) , [1 , 2 * NCPF - 1 , 1]) };
        Sigma_Y(m) = { codistributed(Sigma_Y{m} , dist) };
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    end

    q = 1;
    RMSNorm = cell(1 , cols * (cols + 1) / 2);
    for m = 1 : cols

        for n = m : cols
            RMSNorm(q) = { Sigma_YRef{m} .* Sigma_Y{n} };
            q = q + 1;
        end

    end

        clear  Sigma_Y_U Sigma_Y_V Sigma_YRef_U Sigma_YRef_V 

    %calculating Autocorrelations and Cross-correlation
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %scale factor
    r = 0 : (Length / 2) - 1 ;
    SF1 = (Length / 2) * ((Length / 2) - r).^(-1);
    SF1 = SF1(1 : NCPF);
    r = (Length / 2) + 1 : Length - 1;
    SF2 = (Length / 2) * (r - (Length / 2)).^(-1);
    SF2 = SF2(end - NCPF + 2 : end);
    SF = [SF2 SF1];
    SF = repmat(SF , [NRPF , 1 , WR]);
    SF = codistributed(SF , dist);
    clear r SF1 SF2

    for m = 1 : length(SpectraCell)

        Spec_rho(m) = { Spec_rho{m} / ( NOR * Length ) };
        Spec_rho(m) = { gather(Spec_rho{m}) };
        Spec_rho(m) = { real( ifft( Spec_rho{m} , Length , 2) ) };
        Spec_rho(m) = { [ Spec_rho{m}(: , L0 + NCPF + 1 : end - 1 , :) , Spec_rho{m}(: , 1 : NCPF , :) ] };
        Spec_rho(m) = { codistributed(Spec_rho{m} , dist) };
        Spec_rho(m) = { Spec_rho{m} .* SF };
        Spec_rho(m) = { Spec_rho{m} ./ RMSNorm{m} };
        Spec_rho(m) = { gather(Spec_rho{m}) };

    end
    clear SF RMSNormUU RMSNormVV RMSNormUV


    toc

    disp('saving results')
    tic
    ti = num2str(2 * NCPF - 1);
    tj = num2str(NRPF);
    Rx = repmat(Rx , [NRPF , 1]);
    Ry = repmat(Ry , [1 , 2 * NCPF - 1]);
    Ucell = cell(1 , length(Spec_rho) + 2);
    Ucell(1) = {Rx};
    Ucell(2) = {Ry};
    
    tag1 = num2str(WRInt(1));
    tag2 = num2str(WRInt(2));

    if ispc == 1
        slash = '\';
    else
       slash ='/'; 
    end
    filename = strcat(slash,NameResult,'_Two-point-correlation_',tag1,'-',tag2,'.dat');

    for k = 1 : WR
        zone = num2str(Ry(k , 1));
    %     heading = strcat('VARIABLES = "DX", "DY", "Spec_rho", "Spec_Rvv", "Spec_Ruv", ZONE T = "Yref = ' ,...
    %         zone ,'" I = ' , ti , ' J = ', tj , ' F = POINT');

        heading = strcat('ZONE T = "Yref = ' , zone ,'" I = ' , ti , ' J = ', tj , ' F = POINT');

        for m = 3 : length(Ucell)
            Ucell{m} = Spec_rho{m - 2}(: , : , k);
        end

        data = mixing(2 * NCPF - 1 , NRPF , Ucell);
        data = sortrows(data , [2 , 1]);

        if k ~= 1
            saverappend(folResults,filename,heading,data)
        else
            saver(folResults,filename,heading,data)
        end

    end
    toc
    
    if nn == 1
        NIT = floor(NRPF / WR) + 1;
    end
    
    nn = nn + 1;
    
end
display('Total calculation time')
toc
matlabpool close

